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PACS 75. lO.Nr - Spin-glass and other random models 
PACS 03 . 65 . Ud " Entanglement and quantum nonlocality 
PACS 73 . 43 . Nq - Quantum phase transitions 

Abstract - The entanglement entropy of the random transverse-field Ising model is calculated 
by a numerical implementation of the asymptotically exact strong disorder renormalization group 
method in 2d, 3d and 4d hypercubic lattices for different shapes of the subregion. We find that 
the area law is always satisfied, but there are analytic corrections due to _E-dimensional edges 
(1 < E < d — 2). More interesting is the contribution arising from corners, which is logarithmically 
divergent at the critical point and its prefactor in a given dimension is universal, i.e. independent 
of the form of disorder. 



Introduction. — To study the entanglement proper- 
ties of quantum many body systems is a promising con- 
cept to understand their topological and universal prop- 
erties, in particular in the vicinity of a quantum phase- 
transition point [HE]- Generally the entanglement be- 
tween the subsystem, A and the rest of the system, B, 
in the ground state, is quantified by the von Neumann 
entropy of the reduced density matrix, = TrB|!f')(!f'| as: 
S = — Tr^ (/o_4 log2 Pa)- Generally S scales with the area 
of the interface separating A and B. In some cases, how- 
ever, there are singular corrections to the area law. In one- 
dimensional (d = 1) systems S is logarithmically divergent 
at a quantum critical point [SHS]: 5 = | logg £+cst. Here £ 
is the size of the subsystem and the prefactor is universal, 
c being the central charge of the conformal field theory. 
Recently one considers also generalizations to Rcnyi en- 
tropy and the properties of the entanglement spectrum 

i- 

In higher dimensions our understanding about bipartite 
entanglement is far less complete, the known results are 
almost exclusively about two-dimensional {d = 2) mod- 
els. Considering non-interacting systems, for free bosons 
the area law [7] is found to be satisfied even in gapless 
phases [5] . On the contrary, for gapless free-fermionic sys- 
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tems with short-range lioppings and a finite Fermi surface 
there is a logarithmic factor to the area law 1^. In inter- 
acting d = 2 systems the area law is generally found to 
be satisfied, but in gapless phases and in quantum critical 
points there are additional logarithmic terms, which are 
expected to be universal. This has been demonstrated for 
the d = 2 transverse-field Ising model [TU] and for the an- 
tiferromagnetic Heisenberg model [11]. For the latter the 
logarithmic terms are associated to two sources: i) cor- 
ners on the boundary of the subsystem and ii) non-trivial 
topology in the bulk. There is another class of d = 2 criti- 
cal systems described by d = 2 conformal field theory, the 
prototype being the square lattice quantum dimer model 
[12j . For these models non-perturbative analytical and nu- 
merical results are available and the log-correction to the 
area law is shown to be universal and related to corners 

Besides pure systems there are also investigations about 
the entanglement properties of quantum models in the 
presence of quenched disorder |14| . In d = 1 random sys- 
tems (random antiferromagnetic Heisenberg and XX mod- 
els, random transverse-field Ising model (RTIM), etc.) the 
critical point is controlled by a so called infinite disorder 
fixed point (IDFP) [15], which can be conveniently stud- 
ied by the strong disorder renormalization group (SDRG) 
method [T6l[T7]- Using this approach logarithmic entan- 
glement entropy is found with a universal prefactor |18j . 
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which has been numerically checked by density-matrix 
renormalization jl9j and by free-fermionic methods |20j . 
For ladders of the RTIM the same scaling behavior of the 
entropy is found jH] as in d = 1. In d = 2 the entangle- 
ment entropy of the RTIM has been studied by the SDRG 
method in two papers with conflicting results at the crit- 
ical point. Lin et al |22] has used periodic systems up 
to linear size L = 64 and the numerical results are inter- 
preted in terms of a double-logarithmic factor to the area 
law: 5 ^Inln^. In a subsequent study Yu et al [23] 
has used open systems up to L = 160 and the numerical 
data are fitted with a logarithmic correction to the area 
law: S ^ al + 61n£. This type of choice of the singularity 
is motivated by the similar form of the entropy in = 2 
conformally invariant models [T^, although the logarith- 
mic correction is not attributed to corner effects but to 
percolation of correlated clusters. 

In the present work we revisit the problem of scaling 
of the entanglement entropy in the d = 2 RTIM and use 
an improved numerical algorithm of the SDRG method 
[mUS]. We have studied finite systems up to L = 2048 
and also investigated the entropy of the same model in 
d ~ 3 and d = A. Our goal is to answer the following basic 
questions. 

1. How the criticality of the RTIM is manifested in the 
singular behavior of the entanglement entropy? 

2. What is the physical origin of this singularity, corner 
and/or bulk effects? 

3. Is this singularity universal and independent of the 
form of disorder? 

4. Is it related to the diverging correlation length? 

The rest of the paper is organized as follows. After defi- 
nition of the RTIM we recapitulate the basic steps of the 
SDRG method to calculate the entanglement entropy. Our 
results for c? = 2, 3 and 4 are presented in more details at 
the critical point and afterwards outside the critical point. 
Our Letter is closed by a Discussion and a detailed deriva- 
tion is presented in the Appendix. 

Model and the SDRG method. - The RTIM is 
defined by the Hamiltonian: 

'H = -Y,J.,<r^<r^-Y,h,at , (1) 

{ij) i 

in terms of the erf Pauli-matrices at sites i (or j) of a hy- 
percubic lattice. The nearest neighbor couplings, Jij, and 
the transverse fields, hi, are independent random numbers, 
which are taken from the distributions, p(J) and q{h), re- 
spectively. Following Refs. [24l[25] we have used two dis- 
order distributions, for both the couplings are uniformly 
distributed in [0, 1]. For box-h disorder the distribution of 
the transverse-fields is uniform in [0,/if,], whereas for the 



fixed-h model we have hi = hf, Vz. The quantum con- 
trol parameter is defined as = hi{hb) and 9 = ln{hf), 
respectively. 

The ground state of the RTIM is calculated by an im- 
proved numerical algorithm of the SDRG method [24tf25] . 
During the SDRG [17] the largest local terms in the 
Hamiltonian in Eq.(IT]) are successively eliminated and new 
Hamiltonians are generated through perturbation calcu- 
lation. For the RTIM two types of decimation arc per- 
formed, i) If the decimated term is a strong coupling, say 
Jij, then the two sites, i and j are aggregated to a new 
effective spin cluster, ii) If the decimated term is a strong 
transverse field, say hi, then the site, i is eliminated. After 
decimating all degrees of freedom the ground state of the 
system is found as a collection of independent ferromag- 
netic clusters of various sizes; each cluster being in a GHZ 
state [miin]: ^ (| tt • ■ ■ t) + I U • ■ ■ i))- Each cluster 
contributes by an amount of log2 2 = 1 to the entangle- 
ment entropy if it is shared by the subsystems, otherwise 
the contribution is 0. Thus calculation of the entangle- 
ment entropy for the RTIM is equivalent to a cluster count- 
ing problem. This is illustrated for d = 2 in Fig. [1] (left 
and middle panels). The structure of the GHZ clusters is 
different in the ferromagnetic phase, 9 < 6c, when there 
is a giant cluster and in the paramagnetic phase, 9 > 9c, 
when all clusters have a finite extent. The location of the 
critical point, 9c, has been calculated previously for the 
two random distributions in d = 2 Ref. 1241 and in d = 3 
and 4 in Ref. [251 . 




Fig. 1: (Color online) Renormalized spin clusters (left panel) 
and those connected clusters, which have a contribution to the 
entanglement entropy (middle panel) in a 64 x 64 system with 
a subsystem of 32 x 32 spins (fixed-/i disorder at the critical 
point). (In the left panel clusters with the same mass are repre- 
sented by the same color (greyscale), in the middle panel spins 
in the same cluster are denoted in this way.) Right panel: par- 
tition of the sample into four squares, denoted by 1,2,3 and 
4 and for four slabs, each being composed of two squares as: 
(1, 3); (2, 4); (1, 2) and (3, 4). In this sample and with this sub- 
division there is one "corner" -cluster (see at the left-low corner 
of 1), which contributes to all the four slab subsystems, but 
does not contribute to all the four square subsystems, see text. 
Note that periodic b.c.-s are used and the average value of con- 
nected and corner clusters in this geometry are measured 18.7 
and 0.22, respectively. 

Results at the critical point. — We have calculated 
the entanglement entropy of the RTIM at the critical point 
in finite hypercubic samples of linear size, L with full pe- 
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riodic bomidary conditions (b.c.-s), the largest sizes for 
box-/i (fixed-/i) disorder being 2048 (1024), 128 (64) and 
48 (24) for d = 2, 3 and 4, respectively. (In the latter case 
the clusters are more compact, contain more sites and thus 
the analysis of the entropy is more involved.) We have con- 
sidered d different geometries, in which inO<I?<c?— 1 
directions A extends to the full length of the system, L 
and has periodic b.c.-s, whereas in the other directions 
its length is ^ < L. The three possible geometries for 
d — 3 are illustrated in the inset of Fig. [2] Only in the 
cube geometry with D = there are corners, whereas for 
D = d — 1 in the slab geometry the interface contains no 
edges. For a given random sample and for each geome- 
try we have averaged the entanglement entropy for every 
possible position (and orientation) of A and subsequently 
we have averaged over several samples. The typical value 
of realizations being 40000 but even for the largest sizes 
we had at least 10000 samples. For a given realization the 
extra computational time needed to perform the cluster 
counting problem for the entropy is ©(L^^**"^^), which 
can be speeded up in the slab geometry, see next subsec- 
tion. 

Slab geometry. We start our investigations in the slab 
geometry, where the entanglement entropy of a sample 
averaged over all positions can be written in a simple 
closed form in terms of cluster statistics. Here we an- 
nounce the result, details of its derivation can be found 
in Ref. In this algorithm we consider that axis, say 
the z-axis, which is perpendicular to the surface of the 
slab and we measure the z-coordinate of the points of the 
clusters. For each cluster we arrange the different z values 
as zi < Z2 < ■ ■ ■ < Zk and define the difference between 
consecutive z-values, hi = min[2:i+i ~ Zi,L — (zi+i — Zi)], 
i = 1,2, .. .k - 1; hk = min[zfe - zi,L - {zk - zi)] . Re- 
peating this measurement for all clusters we calculate the 
statistics of the hi differences: n{j) being the number of 
distances with j = h. The position averaged entanglement 
entropy of the sample is then given for £ < L/2: 

41Ui,^) = ^EE"W- (2) 

i=l j=i 

This type of algorithm works in 0{L) times, which is to 
be compared with the performance of the direct cluster 
counting approach: 0{L?). 

In the slab geometry, the area is independent of i and 
any singular contribution to the area law can only be of 
bulk origin. In our study we have fixed L to its largest 
value and calculated the entropy per area, ad-i{L,£) for 
varying £. For 1 ^ £ ^ L wc have found that ad-i{L,£) 
approaches a constant with a correction term: ^ £~'^'^^. 
To illustrate this relation we have calculated the finite dif- 
ference: SS^^tiiLJ) = S^tiLJ + 1) - 41 (i,^) as a 
function of £, which has the behavior: SSlfJ^{L,£) ~ £~'^ 
as shown in Figl2] for d ^ 2,3 and 4. This type of non- 
singular contribution to the entropy in the slab geometry 



can be interpreted in the following way. Due to the finite 
width of the slab only those correlated domains can effec- 
tively contribute to the entropy, which have a finite extent 
^ < f . (Much larger clusters have typically no sites inside 
the slab.) Finite-size corrections are due to clusters with 
£, ^ £, the number of these blobs scales as nti ~ {L/£Y~^ 
and each has the same correction to the entropy, which 
then scales as ^ n^i in agreement with the scaling Ansatz 
and with the numerical data in Figl^] In the SDRG rea- 
soning used in Ref. [22] and which has lead to a lnln£ 
multiplicative correction one assumes the existence of sev- 
eral (^-dependent number of) independent large clusters 
in a ^ ss £ blob, which is in contradiction with the results 
of the present large-scale calculation. 




Fig. 2: (Color online) Finite difference of the entropy in the 
slab geometry, SS^^^^^{L, £), as a function of the width £ for 
d = 2,3 and 4 for the largest sizes, L, for hox-h disorder. (The 
error is smaller than the size of the symbols.) The asymptotic 
forms, ~ £~'^, are indicated by straight lines. Inset: Three 
geometries of the subsystem A studied in d = 3: cube {D = 0), 
column (D — 1) and slab (D = 2). 



Column geometry. In the column geometry, see the 
middle panel in the inset of Fig. [21 there are correc- 
tions to the area law due to edges. Let us consider an 
i?-dimensional edge (1 < E < d — 1) with a total sur- 
face, Je ~ L^, so that its contribution to the entropy is 
given by: aEfs- We have found, that the as prefactors 
have alternating signs: aE/dE-i < and 0^-2 < 0. Here 
using the same reasoning as in the slab geometry the cor- 
rection to the edge contribution per surface is given by: 
oe — 0'e{£) ^ £~^ , which result has been checked numeri- 
cally. Thus we can conclude that the contributions to the 
entropy due to edges are also non-singular and singular 
contributions can only be obtained at corners. 

Cube geometry. In order to check the corner contri- 
butions to the entropy, Scf\£), we study here cube sub- 
systems, as shown in the right panel of the inset of Fig. [21 
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In this case we write >5^ube(^) general case as: 

d-2 

'^iuLW = ad-ifd^i + E '^eIe + s(f\e) . (3) 

E=l 

where the corner contribution has the sign: (— 1)''+^, 
which is opposite to the sign of ai. In d = 2, when 
the subsystem is a square, the second term in Eq.Q is 
missing and we obtain accurate estimates for the cor- 
ner contribution by evahiating the difference: (55^^^ (£) = 
5(2) (^) - 25(2) (^/2) « Si?{i) - 2Si?{i/2). This is pre- 
sented in Fig 13] as a function of In^ for the two types of 
disorder using the largest finite systems. 
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Fig. 3: (Color online) Corrections to the area law in d = 2 
as a function of In^ for fixed-/i L = 1024) and box- 
h (□, L — 2048) disorders. The asymptotic behaviors for 
large £ are indicated by straight lines having the same slopes: 
fo''^' = —0.029. Inset: effective, ^-dependent — fo''^' parameters 
obtained from two-point fits. 

For large £ the data approach a linear logarithmic de- 
pendence, (55(2) (^) _ sS\e) + est ~ -b^^hne + est. 
We have calculated effective, ^-dependent &(2) values from 
two point fits, which are presented in the inset of Fig|31 
From their extrapolation we obtain the estimate, &(2) ~ 
—0.029(1), for both types of disorder, which is to be com- 
pared with 6(2) = -0.019(5) by Yu et al [33] calculated in 
a much smaller system with hax-h disorder. 

In higher dimensions, d > 3, the corner contribution 
represents only a very small fraction of the entanglement 
entropy and thus its estimate through a direct analysis 
of Eq.(|3]) contains rather large errors. We can, how- 
ever, circumvent this problem by considering samples with 
i = L/2, when sir\i) is expressed as appropriate combi- 
nation of the entropies of subsystems with different shapes 
for D = 0, 1, . . . , d — 1. This calculation is presented in 
the Appendix and illustrated in the right panel of FigH] 
for d = 2. Here a given L x L sample contains four square 
subsystems and also four slab subsystems. In the two ge- 
ometries the accumulated boundary between the subsys- 
tems and the environment is the same, thus the difference 



between the accumulated entropies gives the corner con- 

(2) (2) (2) 

tribution: 5cr = 5squarc — S^iah- This contribution is not 
zero, since the so called "corner" clusters, which have no 
sites in one or two non-contacted 90-degree corners (see in 
the right panel of FiglT]), provide different contributions 

to 5square, than tO ^gjg'^^. 

The corner contribution to the entanglement entropy at 
the critical point with £ = L/2 for different values of L 
arc presented in Fig. |3]for d = 2, 3 and 4 and for the two 
disorder distributions. For d = 2 the variation with InL 
is similar to that obtained by the direct analysis in Figl31 
Also for d > 2 asymptotically a logarithmic increase is 
found S^?{£ = L/2) ^ &('^) hv£ + est. and the prefactors 
are estimated by two-point fits. These are shown in 
the inset of FigHl Their extrapolated values are found 
to be disorder independent, thus universal and these are 
listed in the caption of Fig. 2] For d = 2 this coincides 
with the value calculated previously in Figl3l 




Fig. 4: (Color online) Corner contribution to the entanglement 
entropy as a function of In^, for d = 2, 3 and 4 from up to 
down for fixed-/i (-I-) and box-fe (□) disorders. Note that the 
position of of the vertical axis is shifted by 0.1 with d. In 
the inset the effective prefactors of the logarithm are shown 
as calculated by two-point fit. The extrapolated values are 
disorder independent: b'^' = -0.029(1), 6(^) = 0.012(2) and 
fc**^ = -0.006(2). The error of data is smaller than the size of 
the symbols in the main panel and these are smaller than twice 
of the size of the symbols in the inset. The lines through the 
points are guide to the eye. 

The In £ dependence of the corner contributions to the 
entropy can be understood with the example of two-point 
clusters. (Here we remark that during renormalization 
spins are glued together to form new effective spin vari- 
ables and a final spin cluster, which appears in the right 
panel of FiglT] is also the result of the aggregation of two 
effective spins.) It is easy to see, that in d dimension a 
two-point cluster is a "corner" cluster if the two points 
are located in two such hypercubes, which are connected 
by the main diagonal. If the relative coordinates of the 
two-point cluster are Q < Xj < L/2, j = 1, 2, . . . , d (due 
to periodic b.c.-s) then its accumulated contribution to 
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the corner-entropy (obtained by averaging over all pos- 
sible positions) is — 2 J|^^j^(— Xj/i). The probability of 
having a two-site cluster of a length, r, is given by the 
average pair correlation function, Cav('') ~ n'^, where 
Ur ^ r~'^ is the density of non-decimated sites (when the 
typical length between existing effective spins is r). The 
average contribution to the corner-entropy can be esti- 
mated as: Sa\i) ^ ~ dxi . . . /j^ dxd nj=i(~2;j/r^) ^ 
J^^(rd-Vd)/r2d(j^ _ (-l)'*+iln£, which is loga- 
rithmically divergent in any dimension, in agreement with 
the numerical results in Fig. |4l 

Results outside the critical point. We have 
also studied the behavior of the corner-entropy outside 
the critical point and measured Scf\L,d) as a function 
of S = 9 — Oc- In the ordered phase, S < 0, and for 
^ < ^ the giant cluster behaves as a so called global cluster, 
which has 1 contribution to the entropy for all position, 
orientation and shape of the subsystem. As shown in the 
Appendix in odd (even) dimensions, after averaging for 
all positions a global cluster has a contribution 2^"^^ (0) 
to the corner entropy. Approaching the critical point for 
^ > £ these giant clusters have a finite, but ^-dependent 
contribution, that we omit in the following analysis. 

In the upper panels of Fig. \5\ Sif\L, S) is presented 
as a function of S for different finite systems for box-/i 
disorder. For any d the corner-entropy is extremal around 
the critical point and its value outside the critical point is 
well described with the substitution: £—?■£,, with ^ ~ 
being the correlation length. Close to the critical point 
it satisfies the scaling relation: —Scf\L,S) + b'^'^^lnL = 
f{SL^^'^), as illustrated in the lower panels of Fig. [5j Here 
we have used our previous estimates for the correlation 
length critical exponents mUS], v = 1.24, 0.98 and 0.78 
in d = 2, 3 and 4, respectively. 



Conclusion. — We have studied the entanglement en- 
tropy of the RTIM in the vicinity of the quantum critical 
point in dimensions d = 2, 3 and 4 by an efficient nu- 
merical implementation of the SDRG method. Since the 
critical properties of the RTIM are governed by IDFP- 
s [111I151I27] at which the SDRG becomes asymptotically 
exact also our results about the singularities of the entropy 
tend to be exact for large scales. We expect that our finite- 
size results are already in the asymptotic regime, which is 
supported by the fact that the singularity parameters ob- 
tained are disorder independent. We have demonstrated 
that the area law is satisfied for d > 2 and there is a 
singular correction to it in the form of b^'^^ \ni. This cor- 
rection is shown to be attributed to corners, related to the 
diverging correlation length and universal, i.e. disorder 
independent. 

Our investigations can be extended and generalized to 
several directions. First we mention that the extremal be- 
havior of the corner entropy at the critical point makes a 
possibility to detect and define sample dependent critical 
points, a concept which has already been applied in d = 1 
}28j . It is of interest to study the possible singularities of 
the entropy per area and the edge contributions per sur- 
face as a function of 6 around the critical point for d > 2. 
One can also study the entanglement properties of diluted 
transverse-field Ising models for d > 2 having critical prop- 
erties related to classical percolation . Finally, we men- 
tion dynamical aspects of the entanglement entropy after 
a sudden change of the parameters in the Hamiltonian at 
time t = 0. This question has been recently studied [30] 
in d = 1 and an ultraslow increase of the entropy is found: 
S{t) ^ In hit, if the quench is performed to the critical 
point of the system. For d > 2 one expects a similar time- 
dependence of the corner-contribution: Scf\t) ~ Inlni. 
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Fig. 5: (Color online) Upper panels: entropy contribution of 
the corners <si?(L, 5) in d = 2, 3 and 4 for different system sizes 
as the function of the control parameter for hox-h disorder. In 
the lower panels scaling collapse of the central part of <Si?(L, S) 
is shown, see text. 
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Appendix: Corner-contribution to the entropy 

for £ ~ L/2. — Here we show how the corner- 
contribution to the entropy for a given sample can be de- 
duced from the entropies measured in different shapes of 
the subsystems. 

We consider a d-dimensional hypercubic system with 
linear length L with full periodic boundary conditions. 
Inside the hypercubic system we select subsystems of dif- 
ferent shapes, which span the system in D = . . . ,d — \ 
directions, but restricted to length L/2 in the others. (See 
the inset of Fig 15]) The so defined subsystems have hyper- 
faces of dimension A = Z?,Z? + l,...,(i— 1, and the surface 
of a A-dimensional unit is Qd{^) = L^2^~^, while the 
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number of equivalent hyperface units is given by 

71^^ (A) = 2'' 



We measure the entanglement entropy in this system, iSj^'' , 
for all different D, averaged over the L"^^^ possible posi- 
tions and over the (^) orientations of the subsystem. The 
entanglement entropy is written as the sum of the contri- 
butions of the different dimensional hyperfaces: 



A=D 



■'D 



(A), 



(5) 



which is to be compared with Sq"^^ = S^'^^^ in Eq.©. First 
we note that 



to (A); 



its,) 

(1) 



(6) 



where g£)(A)ri,^''(A) and go{A)n'^\A) are the total areas 
of the given hyperface measured in the two shapes. In this 
way we obtain: 



S 



id) 



d-l 



)D \A-DJ 



A=D \a) 



(A) 



(7) 



It is straightforward to check, that this expression can be 
inverted to obtain the entropy contribution in the cube 
geometry: 



d-l 



(-1) 



D-A 



D=A 



2D 



D 



(8) 



As a special case for S^^^^O) = Sc? we obtain for the 
corner contribution: 



5^(0) 



d-l 

E 

D=0 



(9) 



As an application, we calculate the contribution of a global 
cluster to Sq'^\o). A cluster is global by our definition, if 
its entropy contribution is 1 for all positions, orientations 
and shapes (D) of the subsystems, thus 
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(10) 
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